Non-deterministic density classification with diffusive probabilistic cellular automata 
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We present a probabilistic cellular automaton (CA) with two absorbing states which performs 
classification of binary strings in a non-deterministic sense. In a system evolving under this CA rule, 
empty sites become occupied with a probability proportional to the number of occupied sites in the 
neighborhood, while occupied sites become empty with a probability proportional to the number of 
empty sites in the neighborhood. The probability that all sites become eventually occupied is equal 
to the density of occupied sites in the initial string. 

PACS numbers: 05.70.Fh,89.80.+h 



INTRODUCTION 

Cellular automata (CA) and other spatially-extended 
discrete dynamical systems are often used as models of 
complex systems with large number of locally interacting 
components. One of the primary problems encountered 
in constructing such models is the inverse problem: the 
question how to find a local CA rule which would exhibit 
the desired global behavior. 

As a typical representative of the inverse problem, the 
so-called density classification task 1] has been exten- 
sively studied in recent years. The CA performing this 
task should converge to a fixed point of all l's if the ini- 
tial configuration contains more l's than 0's, and to a 
fixed point of all 0's if the converse is true. While it has 
been proved 0] that the two-state rule performing this 
task does not exist, solutions of modified tasks are pos- 
sible if one allows more than one CA rule 0], modifies 
specifications for the final configuration 0, or assumes 
different boundary condition Approximate solutions 
have been studied in the context of genetic algorithms in 
one and two dimensions Q. 

In what follows, we will define a probabilistic CA which 
solves the density classification problem in the stochastic 
sense, meaning that the probability that all sites become 
eventually occupied is equal to the density of occupied 
sites in the initial string. 

We will assume that the dynamics takes place on a one- 
dimensional lattice with periodic boundary conditions. 
Let Si(k) denotes the state of the lattice site i at time k, 
where i G Z, k G N. All operations on spatial indices i 
are assumed to be modulo L, where L is the length of the 
lattice. We will further assume that Si(k) G {0, 1}, and 
we will say that the site i is occupied (empty) at time k 
if 8i(k) = 1 (si(fc) = 0). 

The dynamics of the system can be described as fol- 
lows: empty sites become occupied with a probability 
proportional to the number of occupied sites in the neigh- 
borhood, while occupied sites become empty with a prob- 
ability proportional to the number of empty sites in the 
neighborhood, with all lattice sites updated simultane- 



ously at each time step. To be more precise, let us denote 
by P(si(k + l))\si-i(k),Si(k),ai+i(k)) the probability 
that the site Si(k) with nearest neighbors Si-i(k), s<+i(fc) 
changes its state to Si( k + 1) in a single time step. The 
following set of transition probabilities defines the afore- 
mentioned CA rule: 



P(1|0,0,0) = 
P(1|0,1,0) = l-2p 
P(l|l,0,0)=p 
P(l|l,l,0) = l-p 



P(1|0,0,1) =P 
P(1|0,1,1) = l-p 
P(l|l,0,l) = 2p 
P(l|l,l,l) = l, (1) 



where p G (0, 1/2] (the remaining eight transition proba- 
bilities can be obtained using P(0\a, b, c) = 1— P(l|a, b, c) 
for a,b,c G {0,1}). The probabilistic CA defined by 
JTJ can be defined explicitly if we introduce a set of iid 
random variables {Xi}^ =Q with probability distribution 
P(Xi = 1) = p, P(Xi = 0) = 1 — p, and another set 
{Yi}f =Q with probability distribution P(Y, = 1) = 2p, 
P(lj = 0) = 1 — 2p. Dynamics of the rule |TTJ can then 
be described as 



Si(k + 1) = Xi(l - Si-i)(l - Si)s i+1 
+ {l-Y i ){l-s i - 1 )s i {l-s i+1 ) 
+ (1 - Xi)(l - Si-i)siS i+ i 
+X i s i - 1 (l - Si)(l - s i+ i) 
+YjSi_i(l - Si)s i+ i 

+ (1 - Xi)8i-X8i(l - S l+ i) 

+Si-i8iSi+i. (2) 

To make the above formula easier to read, we omitted the 
time argument, denoting s.j(fc) by s,. After simplification 
and reordering of terms, we obtain 

Si(k + 1) = Si - SiYi + A^Sj-i + Xi.s i+ i (3) 
+ fa-iSi + SjSj+i - 2s. i _is;s i+ i + s i _is i+ i)(y i - 2Xi). 
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DIFFERENCE AND DIFFERENTIAL 
EQUATIONS 

The state of the system at the time k is determined 
by the states of all lattice sites and is described by the 
Boolean random field s(k) = {sj(fc) : i — . . .L}. The 
Boolean field {s(fc) : k — 0, 1,2, . . .} is then a Markov 
stochastic process. Denoting by E B ^ the expectation of 
this Markov process when the initial configuration is s(0) 
we will now define the expected local density of occupied 
sites by Pi(fc) = -E s (o) [ s i(^)]- The expected global density 
will be defined as 



(4) 



i=0 



While both pt(k) and p(k) depend on the initial config- 
uration s(0), we will drop this dependence to simplify 
notation. We will assume that the initial configuration is 
exactly known (non-stochastic), hence p(0) = J2i=o s »(0) 
is the fraction of initially occupied sites. 

Taking expectation value of both sides of © , and tak- 
ing into account that £wo) \Xi — 2X;] = 0, we obtain the 
following difference equation 

Pi(k + 1) = pi{k)+p(p i+ x(k)+pi-i(k) - 2 Pi (k)). (5) 



After summing over all lattice sites this yields 
p(k + l)=p(k), 



(6) 



which means that the expected global density should be 
constant, independently of the value of parameter p and 
independently of the initial configuration s(0). We can 
therefore say that the probabilistic CA defined in (J2J 
is analogous to conservative CA, i.e., deterministic CA 
which preserve the number of occupied sites |jl 1^. Hollll| . 

Note that up to now we have not made any approxima- 
tions, i.e., both JSJ and @) are exact. We can, however, 
consider limiting behaviour of (J5J when the physical dis- 
tance between lattice sites and the size of the time step 
simultaneous ly g o to zero, using a similar procedure as 
described in |l2j. Let x — ei and t — rfc. Now in J^J we 
can replace p(i, k) by p(x, t), p(i ± 1, k) by p(x ± e, t) and 
p(i,k + 1) by p(x,t + r), which results in the following 
equation: 

p(r, f+r) = p(x, t)+p (p(x + e,t) + p(x — e, t) — 2p(x, t)) . 

We will consider diffusive scaling in which time scales 
as a square of the spatial length, meaning that t = e 2 . 
Taking Taylor expansion of the above equation in powers 
of e up to the second order we obtain 



d t p = pd%p, 



(7) 



i.e., the standard diffusion equation. Due to the form 
of |0, in what follows we will refer to the process de- 
fined in J2J) as diffusive probabilistic cellular automaton 
(DPCA). 
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FIG. 1: Fraction of occupied sites N(k)/L as a function 
of time k for two sample trajectories starting from identical 
initial configuration with N(0) = 30, L = 100, p = 0.3. The 
third, almost horizontal line represents average of 1000 such 
trajectories. 



ABSORPTION PROBABILITY 

We will now present some simulation results illustrat- 
ing dynamics of DPCA. Since main features of DPCA 
remain qualitatively the same for all values of the param- 
eter p in the interval (0, 1/2], we have chosen p = 0.25 as 
a representative value to perform all subsequent simula- 
tions. 

Let N(k) — J2i=i s i(k) be the number of occupied sites 
at time k. If we start with N(0) = 0, then N(k) = for 
all k > 0. Similarly, if N(0) = L, then N(k) = L for all 
k > 0. The DPCA has thus two absorbing states, corre- 
sponding to all empty sites (to be referred to as 0) and 
to all occupied sites (to be referred to as 1). If we start 
with < N(0) < L, then the graph of N(k) resembles a 
random walk, as shown in Figure^ Both sample trajec- 
tories shown there eventually end in the absorbing state, 
one of them in 0, another one in 1. This is a general prop- 
erty of the DPCA: regardless of the initial configuration, 
the system sooner or later ends up in one of the two ab- 
sorbing states. Although the time required to reach the 
absorbing state can be large for a given realization of the 
process, the expected value of the number of time steps 
required to reach the absorbing state is finite, as it is 
the case for all finite absorbing Markov chains [13|. Fig- 
ure [21 illustrates this property for L = 100 and the initial 
configuration with 30 occupied sites clustered around the 
center, ie., located at i = 35, 36, 64. All other sites are 
empty. We start with an assembly of 200 of such initial 
configurations, all plotted as vertical lines in which black 
pixels represent occupied sites, while white pixels repre- 
sent empty sites, as in Figure Et- Each of these initial 
configurations evolves according to the DPCA rule, and 
after k — 100 (k — 1000) iterations they are again plot- 
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FIG. 2: Multiple realizations of the CA evolution at (a) 
k = 0, (b) k = 100, (c) k = 1000 and (d) k = 12000. Each 
vertical line corresponds to different realization of the process 
on L = 100 lattice. Black pixels representing occupied sites 
and white pixels empty sites. 200 different realizations of the 
process are shown. 



ted as 200 vertical lines in Figure J2t)- After 12000 
iterations all 200 configurations reach absorbing states, 
as illustrated in Figure EJi. Obviously, some reach the 
state 1, while others 0, yet it turns out that the fraction 
of configurations which ended up in the state 1 is very 
close to 30%, the same as the fraction of occupied sites 
at k = 0. 

To explain this phenomenon, let us define UAr(fc) to be 
the probability that the number of occupied sites at time 
k is N. Since the Markov process {s(k) : k = 0, 1, 2, . . .} 
is finite and absorbing, no matter where the process 
starts, the probability that after k steps it is in an ab- 
sorbing state tends to 1 as A: tends to infinity |l3|. This 



implies that 

lim u N {k) = if N ^ and N ^ L, (8) 

k— >oo 

lim (u L (k) + u Q (k)) = 1. (9) 

k — >oo 

The expected global density, as defined in (@J), is inde- 
pendent of k, hence 

L 

p(0) = L~ x E s (q) [N{k)\ = L- 1 Nu N (k). (10) 

JV=1 

Taking the limit k — > oo of both sides of the above equa- 
tion, and using (jHJ and JJJJ, we obtain 

lim u L (k) =p(0), (11) 

k— »OG 

lim u (k) = 1 - p(0). (12) 

k — >oo 

We have shown that the probability that the DPCA 
reaches the absorbing state 1 is equal to the initial frac- 
tion of occupied sites p{0). The probability that it reaches 
is 1 — p(0). This is in agreement with the behavior ob- 
served in Figure |3 

The above can be viewed as a probabilistic generaliza- 
tion of the density classification process. In the standard 
(deterministic) version of the density classification prob- 
lem we seek a CA rule which would converge to 1 (0) 
if the fraction of occupied sites in the initial string is 
greater (less) than 1/2, i.e., 

lim u L (k) = Q(p(0)), (13) 

lim u (k) = 9(1 - p(Q)), (14) 

k — >OG 

where O(-) is the step function defined as 8(x) = if 
x < 1/2, and G(x) = 1 if x > 1/2. Thus the difference 
between the deterministic and the probabilistic density 
classification introduced here is the replacement of the 
step function in H13I14(I by the identity function in (|1 lt - 

D2J. 

As opposed to detcrministically determined outcome 
in the standard density classification process, in DPCA 
it is just more probable that the system reaches 1 then 
if the fraction of occupied sites in the initial string is 
greater than 1/2, and it is more probable that it reaches 
then 1 if the converse is true. Additionally, DPCA can 
in some sense measure concentration of occupied sites in 
the initial string. If we want to know what is the initial 
density of occupied sites, we need to run DPCA many 
times with the same initial condition until it reaches the 
absorbing state, and observe how frequently it reaches 1. 
This frequency will approximate N(0)/L, with accuracy 
increasing with the number of experiments. 

TIME TO ABSORPTION 

Simulations shown in Figurenindicate that for large k 
the system is typically in a state in which blocks of both 
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FIG. 3: The average time required to reach the absorbing 
state (T) as a function of the initial density p(0) = N(0)/L 
for L = 100, p = 0.3. Each data point (+) represents the aver- 
age taken over 1000 realizations of the process with identical 
initial configurations. The continuous line represents fitted 
parabola (T) = ap(0)(l - p(0)). 



empty and occupied sites are relatively long. We can use 
this observation to obtain the approximate dependence 
of the time required to reach the absorbing state on the 
density of initial configuration. 

If we assume that in a given configuration all occupied 
sites are grouped in a few long continuous blocks, then 
the value oiN(k) cannot change too much in a single time 
step. For simplicity, let us assume that the only allowed 
values of AiV(fc) = N(k + 1) -N(k) are {-2, -1, 0, 1, 2}. 

Since p(k) is time- independent, the probability that 
N(k) increases by a given amount should be equal to the 
probability that it decreases by the same amount in a 
single time step. In the agreement with the above, let us 
define Pj as the probability that AN(k) takes the value 
j, so that pj is non-zero for j E {—2, —1, 0, 1, 2}, where 
P-2 = Pi, P-i = Pi, and 2p 2 + 2p x + p = 1. If T z 
denotes the expected time to reach the state N(k) = 
or N(k) = L starting from N(0) = z, a simple argument 
fl4| yields the difference equation which T z must satisfy 

T z = p 2 T z+2 +piT z+ i +p T z +piT z - 1 +p 2 T z - 2 + l- (15) 

The solution of this equation satisfying boundary condi- 
tions Tq = and Tr, = is given by 



where a — l/(8p 2 + 4pi), meaning that the mean time 
to absorption scales with lattice length as 0(L 2 ). The 
above result would remain valid even if we allowed further 
jumps than ±2 (although the form of the coefficient a 
would be different). 

In order to verify if this result holds for the DPCA, we 
performed a series of numerical experiments, computing 



the average time required to reach the absorbing state 
for 1000 realizations of the DPCA process, for a range of 
initial densities. Results are shown if Figure [3] One can 
clearly see that data points are aligned along a curve of 
parabolic shape, as expected from (TSJ. 



CONCLUSION 

The probabilistic CA introduced in this article solves 
the density classification problem in a non-deterministic 
sense. It is interesting to note that the DPCA conserves 
the average number of occupied sites, similarly as deter- 
ministic rules employed in solutions of related problems 
mentioned in the introduction. Indeed, conservation of 
the number of occupied sites is a necessary condition for 
density classification by CA if one allows modified out- 
put configuration, as recently shown in [l^. This sug- 
gests that a wider class of probabilistic CA conserving 
p(k) might be an useful paradigm in studying how lo- 
cally interacting systems compute global properties, and 
certainly deserves further attention. 
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